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ABSTRACT 



Context. This paper presents a method which can be used to calculate models of the global solar corona from observational data. 
Aims. We present an optimization method for computing nonlinear magnetohydrostatic equilibria in spherical geometry with the aim to obtain 
self-consistent solutions for the coronal magnetic field, the coronal plasma density and plasma pressure using observational data as input. 
Methods. Our code for the self-consistent computation of the coronal magnetic fields and the coronal plasma solves the non-force-free 
magnetohydrostatic equilibria using an optimization method. Previous versions of the code have been used to compute non-linear force- 
free coronal magnetic fields from photospheric measurements in Cartesian and spherical geometry, and magnetostatic-equilibria in Cartesian 
geometry. We test our code with the help of a known analytic 3D equilibrium solution of the magnetohydrostatic equations. The detailed 
comparison between the numerical calculations and the exact equilibrium solutions is made by using magnetic field line plots, plots of density 
and pressure and some of the usual quantitative numerical comparison measures. 

Results. We find that the method reconstructs the equilibrium accurately, with residual forces of the order of the discretisation error of the 
analytic solution. The correlation with the reference solution is better than 99.9% and the magnetic energy is computed accurately with an error 
of < 0.1%. 

Conclusions. We applied the method so far to an analytic test case. We are planning to use this method with real observational data as input as 
soon as possible. 
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1. Introduction 

In the recent past, numerical methods based on optimization 
principles have been used for a number of pro blems associated 



with th e calculation of solar MHD equilibria. IWheatland et al 



d2000h were the first to suggest the use of an optimization 
method to calculate nonlinear force-free fields in the corona 
from photospheric measurements. Since then the optimiza- 
tion method has been extended in various ways, for exam- 
ple by improving cer tain aspects of the original method for 
force-free fields (e.g. Wiegelmannl 2004 ). by introducing ad- 
ditio nal features such as plasma p ressure into the method 
(e.g. Wiegelmann & Neukirch 2006b or by reform ulating the 
method in other geometries (e.g. lWiegermannl2007l) . 

Optimization methods have the advantage of being con- 
ceptua lly straightforward and are reasonably easy to imple- 
ment dWheatland et al.l |2000t Iwiegelmann & Neukirchl [2003 



Inhe ster & Wiegelmannl I20061 see e.g.). At the moment they 
also seem to b e very competitive in terms of computati onal 



also seem to b e very c ompetitive in terms ot computational 
efficiency (e.g ISchrijver et al. \2004. iMetcalfet al.l \200H) . A 



slight disadvantage c ompared to, for e x ample, the Grad- 
Rubin method (e.g. Amari etal. 1997t Wheatland l2006b 
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Inhester & Wiegelmannl 12006c lAmari et alj|2006l) is that they 
still lack the same degree of rigorous mathematical basis exist- 
ing for other methods. 

In the present paper we describe a further extension of the 
optimization method to calculate magnetohydrostatic (MHS) 
equilibria (including pressure and gravity) in spherical geom- 
etry. This is important for calculating global models of the 
corona including information going beyond just the structure 
of the magnetic field. In section [2] we describe the basic equa- 
tions of the optimization method for problem in hand. We then 
give a brief description of the analytical 3D MHS equilibria 
that we use to test the numerical code in section [3] and present 
the test results in section [4] Our conclusions are presented in 
section [5] 

2. Basic equations 

The magnetohydrostatic (MHS) equations are given by 
(VxB)xB-jU Vp-^ pV ( I' = (1) 
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V-B = 0, 



(2) 



where B is the magnetic field, p the plasma pressure, p the 
mass density and *P = — the gravitational potential with 
the gravitational constant G, the solar mass M s and the dis- 
tance from the sun's center r. We do not assume an equa- 
tion of state for the coronal plasma, but leave p and p to be 
independent quantities. To find a magnetic field B, plasma 
pressure p and plasma density p satisfying Eqs. ([]]) and 
©, we follow the spirit of the previous optimization meth- 
ods (e.g.lWheatland et alfeOOOtlWiegelmann & Inhesterll2003 



IWiegelmannll2004t IWiegelmann & Neukirchll2006l) and define 
the functional 



L(B,p,p) 



I 



|(VxB)xB-// Vp-// pV^| 2 



+ IV - B|- 



B 2 



r sin 9 dr d6 d<p. 



(3) 



It is obvious that Eqs. (HJ and (|2]i are satisfied if L — 0. Here 
B is a vector field, but not necessarily a solenoidal magnetic 
field during the iteration. The numerical method is based on an 
iterative scheme to minimize the functional L. To simplify the 
mathematical derivation we define the quantities 

fi a = B 2 [(V x B) x B - p Q Vp - po P VP] (4) 
Q b = B 2 [(V B) B] , (5) 



and rewrite L as 



(6) 



L = j [B 2 Qr a + B 2 Qr h r 2 \ smOdrdOi 

Taking the derivative of L with respect to an iteration parameter 
t, where B, p, p are assumed to depend upon t, we obtain 

1 dL C SB r dp 



X 



WdV- I — -GdS 
v w Js vt 

dp 

— p £l a ■ dS, 

s ot 



where 

F = Vx(£2 a xB)-Q a x(VxB) 



+V(Q b ■ B) - £2 b (V • B) + (n 2 + Q 2 ) B 



and 



G = n x (Q a x B) - n(Q b ■ B). 

Assuming that 

Tt 

dp _ 
— = -v poV ■ Ll a 

ot 

dp GM, 
of r- 



(7) 

(8) 
(9) 

(10) 

(11) 
(12) 



with positive definite parameters p, v and £ and that the bound- 
ary integrals vanish, one can easily see that L is monotonically 



decreasing with t (note that this does not necessarily imply that 
L tends to zero). 

Discretized versions of Eqs. ( fTUb to d!2t . together with ap- 
propriate boundary conditions, form the basis for the numerical 
scheme. The boundary conditions have to be consistent with the 
assumption that the boundary integrals vanish, for example, by 
keeping the magnetic field, pressure and density fixed on the 
boundaries during the iteration. For testing the method in this 
paper we shall take these boundary conditions from the known 
exact solution. For practical applications these would have to 
come from observational data. We remark that due to the intro- 
duction of additional forces the constraints on the consistency 
of the boundary conditions for the magnetic field are some- 
what different from the force-free case. However, the general 
theory of magnetohydrostatic equilibria requires for example 
that the pressure has the same value at both foot points of a 
closed field line under the general conditions assumed in the 
pr esent paper. This is similar to th e Cartesian case discussed 
by I Wiegelmann & Neukirchl (120061) . where the pressure equa- 
tion was forward integrated along field lines using an upwind 
method from one foot point to the other (in the test case de- 
scribed later we make use of the property that the pressure is 
known to be consistent). 

In the numerical implementation based on this method one 
has to choose the product of the time-step At with the three 
numerical parameters p, v and Usually these products have 
to be small enough to achieve convergence and this is ensured 
by an adaptive time-step control. In this paper we have chosen 
the three parameters (multiplied by Af) to have the same values 
on all grid points. Previous experience with applying a simi- 
lar method to force -free magnetic fields in spherical geometry 
( Wiegelmannll2007 ) showed that choosing the same values for 
the entire box can lead to long computing times in the polar re- 
gions due to the distortion of the numerical grid in spherical po- 
lar coordinates towards the poles (note that the poles 6 = and 
9 — n are excluded from the computational domain). This could 
in principle be compensated by allowing for a spatial variation 
of the iteration parameters. 



3. 3D MHS Equilibria 

We use the e xact 3D MHS equ ilibrium in spherical coordinates 
presented by iNeukirch dl994 to test our co de (called case II 
in his p aper). In his paper, Neukirchl (1995) extended earlier 
work bv lBogdan & Low (fl986f on exact 3D MHS equilibria in 
spherical co ordinates. The general method was first found by 
Lowl(ll985l) 'previ ous closely related work can also be found 



iwlQl^QJl) (p revious closely related work can also be round 
Lowl (ll982i and lEowl(ll985ln for Cartesian coordinates and 



developed further in ILowI d 1 99 ll 1 19921 Il993alfbl 12005b . The 



method relies on the presence of an external force (in our case 
gravitation) and the basic assumption that the electric currents 
flow only in surface direction perpendicular to the direction of 
the external force. The additional assumption that the depen- 
dence of the current density on the spatial coordinates has a 
special form involving the magnetic field component along the 
direction of the external force (in our case that is the radial 
component of B leads to a linear equation for that magnetic 
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field component. The plasma density and pressure are deter- 
mi ned from the forc e balance equation. 



by expressions for the pre ssure and the de nsity as given by 
Bogdan & Lowldl986b andlNeukirchl(fl995l). 



Neukirch] d 1995b has extended this method by including a iNeukirchl d 1995b followed [Bogdan & Lowl dl986b and nor 



current density component of the constant-a type (aB with a 
constant), which also allows components of the current density 
in the direction of the external force. It is important to empha- 
size that this does not mean that the total field-aligned current 
density is of the linear force-free type, because the other com- 
ponent of the current de nsity also has a f ield-aligned part. 

The formulation by INeukirchl d 1 995b basically reduces the 
mathematical problem to the solution of an equation similar 
to a Schrodinger equation. A sli ghtly simpler formula t ion of 
the same problem was given by iNeukirch & Rastatterl dl999h 
and some analytical solutions (in Cartesian coordinates) with a 
nonlinear relationship be tween the curren t density and the mag- 
netic field were found by Neukirch d 1997b. A formulation using 
Green's functions was given by lPetrie & Neukirchld2000h. 
Since the magnetic field given in Eqs. (45)-(47) of Neukirch! 



d 1995b contain a number of typogra phical errors, we repeat the 



correct field components here. As in 
the functions 



Neukirchld 1995b we define 



fi(r) = 

Mr) = 
with 



cos q + q sin q 
cos qo + go sin qo 

(3 - q 2 ) cos q + 3q sin q 
(3 - ql) cos q + 3q Q sin q 



q = a(r + a), qo 
We then obtairQ 



a(r + a). 



(13) 
(14) 

(15) 
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These expression tend to the solution presented as case III in 
Bogdan & Low (1986) in the limit a — 0. The magnetic field is 
shown in the top panel of Fig. Q] The solution is completed 



1 Compared to INeukirchl J 1 995t) these expression have been cor- 
rected in the following way: a) a factor 1/2 has been included in the 
second term of B r , b) a factor 1/ sin has been removed from the sec- 
ond term of B^, and c) the sign of the third term of is negative. 



malized the radial coordinate r at a radius of 1.5^0. In the 
present paper, we deviate from this and normalize the radial 
coordinate at a radius of 1 Rq, because we want to impose our 
boundary conditions at the solar surface (= 1 Rq). We then 
carry out our numerical calculation on a box which has an inner 
boundary at 1 Rq and an outer boundary at (1.0 + n/2)RQ « 
2.57 Rq. We use a grid with 20 points in the radial direction, 
40 points in the latitudinal direction and 80 points in the lon- 
gitudinal direction. We exclude the polar regions for numerical 
reasons (see discussion in section |2]i and extend the numerical 
box in only from 11.25° to 168.75°. 

As parameters we choose the parameter a, which deter- 
mines the influence of the non-magnetic forces on the equilib- 
rium, to have the value a = 0.2, and we choose the force-free 
parameter a to have the value a = 0.5. Internally our code 
normalizes the length scale with one solar radius, the magnetic 
field strength with the average radial magnetic field strength on 
the photosphere B ave , the pressure with p = [loplB 2 ^ and the 
mass density with p = fio G M s p/R s B\ 



3 2 

ave 



4. Results 

The numerical code based on the optimization method de- 
scribed in section [2] has been run with initial conditions given 
by a potential field calculated from just the B r component of 
the above exact solution on the boundaries. The initial condi- 
tions for density and pressure are as in a stratified atmosphere. 
This configuration balances the radial pressure gradient and the 
gravity force. Consequently the initial plasma is structured only 
in the radial direction and p and p are invariant in and <p. We 
use the exact solution to prescribe the boundary conditions for 
B, p,p. Our code solves the magnetohydrostatic equations (HJ- 
d2J in the computational box with respect to these boundary 
conditions. To evaluate the performance of our code we com- 
pare the result with the exact solution. 

A visual impression of the exact magnetic field is given in 
the top panel of Fig.Q] For comparison the potential field used 
as initial condition for the iteration is shown in the middle panel 
of Fig. Q] One notices in particular that the potential field has 
a different connectivity than the analytical MHS solution. The 
numerical method will thus have to change the magnetic field 
connectivity during the iteration process. 

In Fig. [2] we show r-0 cuts of pressure (left) and density 
(right) for the angle <p — 0. The top panels show these quantities 
for the exact solution, the middle panels for the initial condition 
and the bottom panels show the result of the numerical calcula- 
tion. Visual inspection of these figures gives the impression that 
the method does achieve a good, but not perfect agreement with 
the exact solutions. Especially the regions close to the bound- 
aries still show noticeable differences, which may be caused by 
the general problems with convergence closer to the poles due 
to the deformation of the numerical grid using spherical coor- 
dinates. A more sophisticated numerical grid, as discussed in 
section|5] will probably help to overcome these problems. 
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To check the quality of the reconstruction in a more quan- 
titative way we use a number of methods. The results are pre- 
sented in table Q] First, we evaluate how well the force balance 
condition and the solenoidal condition are fulfilled. To do this 
we evaluate 

- the functional L as defined in ([3), 

- the functional U = f l< VxB >xW^° ™ P f ^ sin g dr d6d<p, 
telling us how well the force balance condition is satisfied, 

- the functional Lq = f |V ■ B| 2 r 2 sin9drd9d(/>, which tells 
us how well the solenoidal condition is satisfied. 

The evolution of the functionals L, L\, Lo during the numeri- 
cal computation is shown in figure [3] in a double logarithmic 
plot. The functionals L and L\ decrease rapidly by about five 
orders of magnitude. The functional L%, which represents the 
solenoidal condition decreases by about three orders of magni- 
tude. The entries in the table give the values of the functionals 
at the end of the numerical calculation. 

(2006) 



Following 



Amari et al 



and 



Wiegelmann & Neukirch ( 20061) we also provide the in- 
finity norms 



- || V ■ B 

- lljxB 



Vp - pvy 



which are defined as the supremum of the divergence and the 
total force density. For comparison, we also provide these val- 
ues for the exact analytic solution if used in a discretized form. 
This allows us to estimate the discretisation error introduced 
by calculating the solution on a finite numerical grid (Here 
20, 40, 80 grid-points in the r, 9, (p direction, respectively). 

We compare the reconstructed equilibrium directly with the 
known analytic eq uilibrium by a number of quantitative mea- 
sures as defined bv lSchriiver et al. I J2006L see section 4 of that 
paper) . The figures quantify the difference between two dis- 
cretized vector fields B (known analytic field evaluated on the 
numerical grid) and b (reconstructed field). These measures 
are: 

- the vector correlation C vec = 2/ Bi V (£ ; |Bi| 2 |bi| 2 ) 1/2 , 

- the Cauchy-Schwarz inequality Ccs — jj 2; |B-'|'|b'| ' wnere N 
is the number of vectors in the field; 

- the normalized vector error £n = |bi - Bj|/ |Bi|, 

- the mean vector error Em = 2; 

- the magnetic energy of the reconstructed field divided by 



the energy of the analytical field e 



S, |b t |- 



Z, IBiP " 

The quality of the reconstruction of the pressure p and the den- 
sity p is quantitatively assessed by correlating the analytic and 
the reconstructed solutions using the linear Pearson correlation 
coefficients (called correlation p and correlation p in table [TJ. 
Finally, we provide the number of iteration steps and comput- 
ing time for a single processor run on a common workstation. 



Table 1. The table provides several figures of merit which 
can be used to assess the quality of the reconstructed solution. 
We compute all figures for the complete computational do- 
main. The analytical reference field was specified in the cones 
11.25° < 9 < 168.75°. For the calculation presented here we 
used the following iteration parameters: p — 1, v = 0.01, £ 
I 





Ref. 


Potential 


Reconstruction 


L 


0.003 


0.004 


0.002 


u 


0.001 


0.002 


0.001 


u. 


0.002 


0.002 


0.001 


II V • B |U 


0.448 


0.582 


0.448 


||JxB-Vp-pVP|L 


0.137 


0.186 


0.143 


r 


1 


0.946 


0.9997 


Ccs 


1 


0.810 


0.997 


E N 





0.378 


0.021 


Em 





0.524 


0.042 


£ 


1 


0.951 


1.0008 


Correlation p 


1 


0.975 


0.9998 


Correlation p 


1 


0.928 


0.9990 


No. of Steps 






214220 


Computing time 






YlKllmin 



d Wheatland et al. I l2000l) to global magnetohydrostatic equilib- 
ria including the pressure force and the gravitational force in 
spherical geometry. The proposed generalization of the opti- 
mization method leads to two additional equations for the pres- 
sure and the density that have to be solved simultaneously to the 
magnetic field equation. Boundary conditions for the magnetic 
field, the pressure and the density are necessary to complete the 
problem. 

We have implemented a numerical code based on the pro- 
posed method and have tested the code usin g a known three - 
dimensional magnetohydrostatic equilibrium ( Neukirchll995 ). 
The numerical calculation is started from a potential field with 
the same radial magnetic field component as the analytical 
equilibrium on the boundary. The initial pressure and density 
distribution are a spherically symmetric stratified atmosphere 
in hydrostatic balance. Both visual inspection of the results as 
well as a quantitative analysis using various diagnostic mea- 
sures indicate that the method works well and converges to the 
analytic equilibrium. For the presented tests we used a low spa- 
tial resolution and got a relatively long computing time (about 
200 000 iteration steps) until convergence. In exp eriments with 
the fo rce-free version of our spherical code (see IWiegelmann 



20071) we found that the computing time scales with N re- 



garding the number of grid points N in one spatial direction. 
This is somewhat slower as the theoretical estimate of A^ 5 



for a c artesian optimization code obtained by I Wheatland et al 
J2000l) . The spherical magnetohydrostatic code is significantly 
slower than the cartesian force-free code for two reasons. 



5. Conclusions and Outlook 

We have extended the optimization method originally pro- 
posed for the reconstruction of force-free magnetic fields 



1. The convergence of the numerical grid towards the poles 
requires a sufficiently small time-step. 

2. The plasma f3 might vary strongly in the entire region and 
in particular low-/3-regions require very small time-steps to 
compute the magnetic field and plasma simultaneously, be- 
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Fig. 4. Outlook: How can the tool become applied to data? The basic idea is to compute first a nonlinear force-free field from 
observed vector magnetograms and then model the plasma along the loops. Our newly developed magneto-hydrostatic code uses 
the resulting magnetic field and plasma configuration as input for the computation of a self-consistent equilibrium. 



cause small changes in the Lorentz force can result in con- 
siderably large changes in the low f3 plasma. 

Point 1.) can be addressed by using a more sophisti- 
cated numerical g r id, e.g . the Yin- Yang grid developed by 
Kagevama & Sato) d2004), which has been applied in geo- 



physics (see e.g. Yoshida & Kageyamal2004 ). This overset grid 
contains two complementary grids which lead to an almost uni- 
formly spaced spherical grid. An additional advantage of the 
Yin- Yang grid is that it is suitable for massive parallelization. 
The Yin- Yang grid has been applied for geophysical simula- 
tions on the Earth simulator super-computer in Yokohama. To 
speed up the 2. point one might compute first the magnetic 
field alone as a nonlinear force-free field (which is a reason- 
able approximation in low-/3-regions) and switch on the self- 
consistent plasma iteration only after for a fine-tuning. A multi- 
scale approach, as recently implement ed to speed up our fo rce- 
free cartesian optimization code (see iMetcalf et all 120071 for 
details) is also an option worth trying for the spherical imple- 
mentation of our force-free and magnetohydrostatic codes. We 
are confident, that the above mentioned potential for improve- 



ments together with a massive parallelization will allow us to 
apply our newly developed method to real data with a reason- 
able grid resolution. 

In figure [4] we outline a scheme on how the code might 
be applied to data. As boundary conditions on pressure and 
density are not directly measured, we propose the follow- 
ing approach. The basic idea is to compute first a force-free 
magnetic field and then model the plasma along the mag- 
netic loops, e.g., by the use of scaling laws and optionally 
with the help of a tomography cod e. Such an approach has 
been used by ISchrijver et all (I2004h by using a global poten- 



tial magnetic field and specifying free scaling law parameters 
(e.g., heating function) by comparing artificial plasma images 
(created by line-of-sight integration from the model plasma) 
with X-ray and EUV observations. We propose to generalize 
this approach by using a nonlinear force-free magnetic field 
model and compare the model plasma with observations from 
two viewpoints as provided by the STEREO-mission. Optional 
STEREO-images can be used additionally to approximate the 
coronal density distribution by a tomographic inversion. As a 
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consequence of this step (modelling the plasma along a mag- 
netic loop) the plasma pressure is consistent along the loops. 
Different values for the pressure on different field-lines will vi- 
olate the force-free condition for the magnetic field, however, 
and the configuration is not exactly in a magnetohydrostatic 
equilibrium. Finite pressure gradients have to be compensated 
by a Lorentz force. This computation can be done with the help 
of the program described in this paper. We propose to use the 
force-free magnetic field configuration and the model plasma 
as initial state for our newly developed magnetohydrostatic 
code to compute a self-consistent MHS -equilibrium. For a low 
P plasma the back-reaction of the plasma onto the magnetic 
field will be small, for higher values of /3 (as found e.g., in hel- 
met streamers) the magnetic field might change significantly. 
As a result of this approach one has reconstructed the 3D coro- 
nal magnetic field and plasma configuration self-consistently 
within the magnetohydrostatic approach and the model is con- 
sistent with measured photospheric vector magnetograms and 
observed coronal images from different viewpoints as well. 
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Wiegelmann et al.: Spherical optimization code 
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Fig. 2. Plots of the spatial variation of the pressure p (left col- 
umn) and density p (right column) in logarithmic scaling in the 
r-6*-plane for = 0. The top, center and bottom panels corre- 
spond to the exact solution, the spherically symmetric stratified 
atmosphere used as initial state for the iteration (also identical 
to the background atmosphere used in the analytical solution) 
and the result of the numerical calculation. The main features 
of the analytical solutions are clearly recovered by the numeri- 
cal method, but some small differences are still present. 
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Fig. 1. Magnetic field lines for the exact analytical solution, a 
potential field with the same B r on the photospheric boundary 
(r = Rq) and the magnetic field obtained by the optimization 
method. The color coding corresponds to the value of B r on 
the photosphere (yellow:positive, blue: negative) and the disk 
center corresponds to <p — 180°. The potential field is used as 
the initial field for the numerical calculation and clearly has a 
different connectivity from the exact analytical solution. The 
reconstructed solution matches the analytical solution down to 
plotting precision, except for one equatorial field line. 
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Fig. 3. Evolution of the functionals L (solid line), L\ (dotted 
line) and L2 (dashed line; for definitions see text) during the 
iteration process. All three functionals decrease by several or- 
ders of magnitude during the iteration. 



